Localized modes in arrays of boson-fermion mixtures. 
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It is shown that the mean-field description of a boson-fermion mixture with a dominating fermionic 
component, loaded in a one-dimensional optical lattice, is reduced to the nonlinear Schrodinger 
equation with a periodic potential and periodic nonlinearity. In such a system there exist localized 
modes having peculiar properties. In particular, for some regions of parameters there exists a lower 
bound for a number of bosons necessary for creation of a mode, while for other domains small 
amplitude gap solitons are not available in vicinity of either of the gap edges. We found that the 
lowest branch of the symmetric solution either does not exist or exist only for a restricted range of 
energies in a gap, unlike in pure bosonic condensates. The simplest bifurcations of the modes are 
shown and stability of the modes is verified numerically. 
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I. INTRODUCTION 

Localized modes constitute an intrinsic feature of non- 
linear systems with periodically varying parameters 0. 
They also represent a signature of other fundamental 
physical phenomena like instabilities and phase transi- 
tions, in particular transition between superconducting 
and insulating states in condensates of atomic vapors. 
That is why during the last few years nonlinear modes 
in Bose-Einstein condensates (BECs) embedded in opti- 
cal lattices, attracted a great deal of attention (see e.g. 
0, 0for review) . Such modes were found to exist in sin- 
gle jj] and multicomponent _5J condensates, as well as in 
atomic-molecular condensates The properties of lo- 
calized modes are governed by the lattice parameters and 
independently by the number of atoms which determines 
position of the chemical potential in a gap of the lattice 
spectrum. 

In the present paper we report the existence of local- 
ized modes in boson-fermion mixtures with large number 
of spin-polarized, and thus noninteracting, fermions. The 
main features of these systems steam from the fact that 
fermionic component is linear and at the same time mod- 
ifies linear and nonlinear properties of the effective me- 
dia for bosons. The role of fermions, when their number 
significantly exceeds the number of bosons, is of crucial 
importance for stability of mixtures 0, 13 > as well as for 
possibility of existence of quasi-one-dimensional solitary 
waves in spatially homogeneous traps and gap soli- 
tons in the presence of the optical lattice 10]. The gap 
solitons, considered in Ref. jlOj were however restricted 
to small amplitudes and to the case where the fermionic 
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component does not result in the spatial variation of the 
nonlinearity. At the same time it was also shown that 
when the Fermi energy is of order of the amplitude of the 
lattice potential, it becomes strongly dependent on the 
spatial coordinate. If in such situation the boson-fermion 
interaction is not negligible compared with boson-boson 
interaction, then in the mean-filed approximation, the 
fermionic component significantly changes not only the 
linear potential but also the effective two-body interac- 
tions among bosons. 

Departing from the known results and taking into ac- 
count that the fermionic distribution itself is determined 
by the trap potential, one can predict that intrinsic lo- 
calized modes in boson-fermion mixtures can exist, and 
that they can possess peculiar properties originated by 
fermions modifying the effective lattice potential and, 
what is most importantly, introducing a nonlinear lattice 
for bosons, i.e. periodic modulation of the nonlinearity 
governing boson-boson interactions. Description of such 
modes and of their properties is our main goal. 

More specifically, in Sec. II we derive the nonlinear 
Schrodinger equation with periodic potential and period- 
ically varying nonlinearity, as a model describing quasi- 
one-dimensional boson-fermion mixture in an optical lat- 
tice. The diversity of localized modes and their proper- 
ties are described in Sec. III. The outcomes of the work 
are summarized in Conclusion. 



II. EVOLUTION EQUATIONS 
A. Mean-field approximation 

We consider low-density bosons and spin-polarized 
fermions at large density embedded in an optical lattice. 
At the zero temperature, the dynamics of fermions in the 
vicinity of the Fermi surface, is described in the hydrody- 
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namic approximation Boson- fermion interactions, 

which are relatively weak due to small density of bosons, 
can be accounted as corrections to the Gross-Pitaevskii 
(GP) equation for bosons. The respective mean-field 
equations were derived in 0- In the present subsection 
we recover the main result of , by means of alternative 
approach, based on the kinetic theory 'TT'|. 

To this end we start with the Hamiltonian of the boson- 
fermion interactions Hint = gbf J rfr\['^^\l/$^$, where $ 
and are the annihilation field operators for bosons 
and fermions, correspondingly, gbf — 2Trh?abf /m, m — 
rubrnf /{nib + '^ij), abf is a boson- fermion scattering 
length, and mbj are atomic masses (subindexes b and 
/ stand for bosons and fermions, respectively). Next 
we introduce the expectation values: the order param- 
eter of bosons ^ = (^) and the averaged density of 
fermions p — (<ji>t$). Then after the standard approx- 
imation \1/ « the equation governing the dynamics of 
bosons acquire the form 



in-— = — - — 

ot 2mb 



where gbb = Airh^abb/mb and abb is the scattering length 
of boson-boson interactions. Hereafter Vbj{r) stand for 
the trap potentials of the two components. 

It follows from the definition of Hint and from Eq. 
(0) that p{r,t) and |^(r,i)p can be interpreted as ex- 
ternal time-dependent potentials applied to bosons and 
to fermions, respectively. This in particular means that 
F = — Vj^'p is the average force which bosons exert on 
a fermion. Hence, the kinetic equation for the distri- 
bution function of the fermions n(r, p, t) has the form 

(Vp = d/dp) 



VpnVe - 



FVpn = . 



(2) 



Here e = £(r, p, t) is the energy of a fermion and link 
between the distribution function and average density of 
fermions is given by 



y n(r,p,t) 



dp 



Following the standard procedure (see e.g. 0), con- 
sidering zero temperature, we represent the distribution 
function in a form of an unperturbed part no(e) depen- 
dent on the unperturbed Fermi distribution and consid- 
ered as a function of the particle energy e, and its exci- 
tation Sn{r,p,t): n = no{e) + 5n(r,p,t). For the spin- 
polarized and thus noninteracting fermions in an external 
potential, we have e = p^/(2m) -I- V/(r). Now we define 
the averaged momentum P{r,t) and the current j(r, i) 
respectively by the formulas 



dp 



(27r?i)3' 



where v = Vp£ = p/ruf is the velocity of a fermion. 
Taking into account that Ve = VV/ and rewriting 

dp 



Sn 



{2TTh)- 



-,(3) 



we obtain in the leading order (a, f3 = x,y, z) 



dt 



dxo 



dp 



+9bf 



{2TThy 

dp 



13 
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dxp 



Sn 



(4) 



In the right hand side of this equation we have accounted 
only the leading terms. In the case at hand j = P/nif 
what allows us, after differentiating with respect to time 
and retaining the terms of the leading order, to rewrite 
the continuity equation 



A^f -I- H(r)1' (1) as follows 



dp 
dt 
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Next we approach 



dp 



3 dx/s 



PfPi, 



where pf is a momentum of the Fermi surface: p\ = 
?i^(67r^/9o)^'''^jand use the stationary Thomas- Fermi dis- 
tribution nHi) 

Po-[2TO/(i?^-y/(r))]^/V(6^'ft'). (5) 
Combining these approximations with we obtain 



9V 

a<2 



= V 



(6) 



Eqs. and © make up a system governing the 

mean-field dynamics of boson-fermion mixture (from the 
direct averaging of the equations for the field operators, 
this equation was earlier obtained in Ref. (3|). 

In the present paper we are interested in specific so- 
lutions, representing localized excitations whose spatial 
dimension is of order of the lattice constant (contrary 
to gap solitons, considered in which extend to many 
lattice periods). This is the case where the stationary ap- 
proximation dpi/dt = is valid. Then © immediately 
gives the relation 



Pi = 



_ ^gbfmf 1/3 



(6^2)2/3^2 



Po' m 



(7) 
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Substituting this formula in the expression for p and sub- 
sequently in (Q, we arrive at 

ih'j- = -—/:^^ + v{v)^ + g{vmH (8) 

ot Zmii 



where 

V{v) = V6(r) + .gb/po(r) 
is the effective lattice potential and 



= gbb 



Q 2 1/3 

(67r2)2/3;j2 



(9) 



(10) 



is the effective nonlinearity. 



B. One-dimensional limit 

Consider now a one-dimensional (ID) lattice along the 
X-axis and tight confinement in the transverse direction: 

Vbj = + Ubjinx). 

Here uj^j are the linear oscillator frequencies in the trans- 
verse direction, = (y, z), k = Tr/d, d is the lattice con- 
stant, and Ubj{nx) are 7r-periodic functions: Ubjinx) — 
Ubjinx + tt). Pauli's exclusion principle leads to sig- 
nificant transverse extension of the distribution jOo(r). 
However, as it is shown in strong confinement of the 
bosonic component and weakness of the two-body inter- 
actions allow transition to the simplified ID evolution 
equation. Thus we require a d, where a — \J h/mhOJb 
is the transverse linear oscillator length of bosons. The 
self-consistent way of the derivation of the ID reduction is 
based on the multiple-scale expansion, which in the lead- 
ing order yields ^'(r, t) = C(r^)^'(a;, t) exp {—iustt), where 
= a~^7r~^/^ exp {—r'^/ (2a^)) describes the linear trans- 
verse distribution and ^'(a;, t) is a smooth envelope. Since 
the respective procedure was described in a numerous 
works (see e.g. for fcrmion-boson mixtures) we skip 
the details. The result is a nonlinear Schrodinger (NLS) 
equation with a periodic potential and periodic nonlin- 
earity, which in dimensionless variables X = kx, T = 
ij{X,T) = {KNb)-^/^^{x,t), and g = po/n^, 

reads 



5^ 



-u{x)i^ + g{xM''^. 



(11) 



The periodic coefficients are obtained from ^ and H1(J|I . 
They are explicitly determined by the stationary distri- 
bution of fermions g{X): 



U{X)=Uo{X)+UigiX), 

g{x) = go-gig'^Hx), 



(12) 
(13) 



gi = 2(6/7r)^/^ {abf/a)'^(mfmb/m'^)Nb. The both func- 
tions U{X) and g{X) are 7r-periodic: U{X) = U{X + it) 
and g{X) = g{X + -k). To make all values to be of the 
unity order, above we introduced Nb, which determines 
the order of magnitude of the total number of bosons. 



III. LOCALIZED MODES 

A. Modulational instability and bifurcation of 
localized modes from the continuum spectrum 

We start the analysis of by the study of the stabil- 
ity of the linear Bloch states with respect to smooth mod- 
ulations of their amplitude. This phenomenon is referred 
to as modulational instability (or alternatively as dynam- 
ical instability) and provides the information about small 
amplitude gap solitons. We designate by s'^a'^ and by 
Lp"a\x) the energy and the Bloch states of lower ("ct" 
stands for "— ") and upper ("cr" stands for "-I-") edges of 
the a's band (a = 1, 2...) of the spectrum of the operator 
C = -d'^/dX'^ +U{X): Llp^S^ = Then the in- 

terval refers to the a's gap and a = corre- 

sponds to the semi-infinite gap. Next, we compute the in- 
verse effective mass {m\^\^ = {l/2)d^Ei"'^ /de and the 

nonlinearity coefficient Xa'' = g{X)[ip^a\xy\^dX . 
Finally, we formulate the condition of the modulational 
instability, which is also the condition necessary for exis- 
tence of small amplitude gap solitons, as 



(14) 



UoiX) = Ub{X)/ER, Er = h^K^/2mb is the boson re- 
coil energy, Ui = AiTKabfmb/rn, 5o = 4^abbNb / ko^ , and 



(see e.g. 0,13 for more details). 

The main feature of the problem at hand, making it 
very different from the standard and well studied case of 
the nonlinear Schrodinger equation with a periodic po- 
tential, is that the obtained criterion (|14|l strongly de- 
pends on the fermionic distribution responsible for the 
spatial dependence of g{X) [see Eqs. (|12ll . (|13|) ] not only 
through the effective mass , but also through the ef- 
fective nonlinearity Xt^ ^ ■ In its turn, the spatial distribu- 
tion g{X) strongly depends on the transverse frequency 
LUf (particular numerical examples can be found in 10]). 
On the other hand presence of the linear potential makes 
the problem very different from the situation studied re- 
cently in Refs. [l3.[l3l| where the spatially periodic non- 
linearity was considered in absence of the linear potential. 
Namely, in a linearly homogeneous, U{X) = 0, nonlin- 
ear lattice one has ip''a\x) = 1 and thus Xa^ — Gav, 
where gav is the spatial average of g{X). Thus Xcf ^ = 
when gav — and no small amplitude solitons can exist 
near gap edges 12]. This leads to existence of a lower 
bound for the number of particles necessary for excita- 
tion of a localized mode. In the presence of a linear lat- 
tice, U{X) 7^ 0, generally speaking XcT^ 7^ gav and even 
signs of Xti^ and gav niay be different. Then, if H14|) is 
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FIG. 1: (Color online) (a) The two lowest bands of the ef- 
fective band structure for bosons and (b) dependence Xq'' 
on Kahh for a mixture of N j = 10* of K*'' atoms per one 
ceU and = 5000 of **'^Rb atoms, at na^j f» 0.0314. The 
lattices are Wb = -10cos(2X) and Uj = -21.75 cos(2JQ 
is taken into account that Uf /Ub ~ rnt/mf = 2.175 
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1.4359 X 10 indicate values of natb where the nonlinearity 
coefficients become zero. 



satisfied, the lower boundary of the number of atoms is 
removed even for Mi^^Gav > 0. 

In Fig^ b we show the nonhnearity coefficients cor- 
responding to the semi-infinite, Xi~\ ^^^d first lowest, 

Xi^'' and X2 ^ gE^PS for the ®^Rb - ^"K mixture, in the 
region where they change the sign. It follows form the 
figure, that small amplitude gap solitons can be created 
neither in the semi- infinite gap for KUbb > 7i \ nor 
in the vicinity of the upper edge of the first band for 
Kflbf, < 7{^-' , nor in the vicinity of the lowest edge of the 
second band at nabb > 72 • An unusual situation occurs 

for 72 < KUbb < 7^"''^: gap solitons do not exist in the 
vicinity of either of the edges of the first gap. 



B. Localized modes in the semi-infinite gap 



We look for a stationary solution of Eq. (|ll|l in the 



form 



-i£T 



4>{X), where 4'{X) is real and (/) ^ as 



X 00. Starting with the semi-infinite gap, i.e. with 
£ < £[~\ we observe, that if Gm = min^c G{X) < 0, then 
in the limit £ —f —00, there exists a soliton, which is 



strongly localized about Gm- More specifically this hap- 
pens when that soliton width, which is estimated as 1 /\/£ 
is much less than the lattice period, i.e. £ 3> 1. Now the 
periodic potential can be viewed as a perturbation and 
the respective solution can be approximated by 



0s (^) 



y/2\£~Urn\/\Gn 



cosh(X-\/|5 — U„ 



(15) 



where Km — minx U{X). 

This conclusion is confirmed by numerical study re- 
ported in Fig. 121 When Xi < Oj the normalized number 
of bosons, defined a.s N — J (j)-^{x)dx is a decreasing func- 
tion of £ (the line Xi ^ = —5.37 in the upper panels). It 
approaches zero as the energy detuning decreases accord- 
ing to the law N cx \J £[^'' — £. This is the standard sit- 
uation 0. If, however Xi ^ > 0, N achieves its minimum 
Nm ~ l.lOSNb in the vicinity of the gap edge and starts 
to grow with the decrease of the energy detuning (the 

line Xi ' = 0.41 in the upper panels). Thus there exists 
a nonzero minimal number of bosons necessary for cre- 
ation of a localized mode (in linearly homogeneous NLS 
equation with a periodic nonlinearity this effect was ob- 
served in E2)i what corroborates with the above conclu- 
sion about impossibility of existence of small-amplitude 
gap solitons at Al[ ^Xi ■* > 0- This phenomenon mani- 
fests itself in shapes of the localized modes: c.f. the pro- 
files in panels A and B, corresponding to the same energy. 
The solutions in Panels C and D display approach to the 
limit 0s (X) obtained above. From the upper panels of 
FigElone can see that for equal energies the number of 



Xi 



> exceeds one at 



the bosons in localized modes at 

x[~^ < 0. 

By direct numerical solution of Eq. (|ll|l we have ver- 
ified that the licalized modes displayed in the panels A, 
B, C and D of Fig.[21are stable. More specifically we per- 
turbed the mode profiles by the factor 1 -I- 0.1 cos (21 X) 
and computed the dynamics until T = 500, observing the 
oscillations of the soliton amplitudes of order of 10% of 
their averaged values. 



C. Gap solitons in the first gap 

Considering the localized modes of the first gap (see 
Fig. OJ, the first intriguing property which one can ob- 
serve at Kttbb > 7i~''^ , i-e. for the parameters giving 
Xi^"* > and xi'^ > 0, is a zig-zag type dependence of 
the boson number on the energy (Fig. |3^). The depen- 
dence of this type can be viewed as a set of successive bi- 
furcations of the branches in the points £i^'^ (j = 1, 2, ...), 
the first two indicated in Fig. Moreover there ex- 
ists a critical energy, coinciding with the first bifurcation 
point £i^\ such that no lowest-branch gap solitons exist 



at < £ < £ 



(-) 



(according to our numerical results 
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FIG. 2: The number of bosons, N/Nb, of the locahzed modes vs energy £ in the semi-infinite gap —oo < £ < £{ (upper panels; 
f}^' coincides with the right boundary of the panel) for the ^^Rb-^^K mixture with the same parameters in Fig. except 
natt ~ 1.4445 x 10~^ (the curves Xi ~ —5.37) and ttatb = 1.4558 X 10"^ (the curves Xi ~ 0.41). The lower panels show 
solitonic shapes corresponding to the points A, B, C, and D (grey and white colors correspond to half-periods with U[X) < 
and U[X) > 0). Only the lowest branch of the solutions is shown. 



similar statement applies also for upper branches of the 
solutions). Motion along the curve N{8) is accompanied 
by the redistribution of atoms among the potential min- 
ima as it is illustrated in Fig. |3|:,d. The unlimited grows 
of N can be understood as follows. When the energy de- 
tuning toward the gap grows, the amplitude of the mode 
grows and the width decreases. Such a mode is localized 
in the vicinity Qm < 0. Increase of the number of parti- 
cles leads also to occupation of the regions where not only 
U{X) > 0, but also G{X) > 0, i.e. where the nonhnearity 
enhances the repelling of the atoms, thus allowing storage 
of a larger number of atoms, compared with nonlinear ly 
homogeneous models. The modes A and B in Fig. I^t- 
were found to be dynamically stable. As in the previous 
case we integrated the Ea. (|ll() with the initial condition 
corresponding to the modes A and B perturbed by the 
factor 1 -I- 0.1 cos (21 X). We observed the oscillations of 
the mode amplitudes of order of ±5% until T — 500. At 
the same time the numerical simulations evidenced, that 
modes C and D are not stable. Their dynamics is accom- 
panied by the travelling of the atoms among minima of 
the lattice. 

Another situation which does not exist in nonlin- 
ear ly homogeneous structures, corresponds to the inter- 
val 7^^^ < Katb < 72 ^ where Xi^"* < and X2~^ > (see 
Fig. Oh) • Here we concentrate on antisymmetric excita- 
tions. Similarly to the case of the semi-infinite gap, there 
exists a minimal number of bosons, which is necessary for 
creation of a localized mode: the lowest branch depicted 
in Fig. does not reach zero. 

By direct simulation of Eq. (|ll|l we found that so- 
lutions corresponding to the part of the branch, where 
dN/d£ < [e.g. modes E and F shown in Fig. Oj], 



are dynamically stable, while solutions corresponding to 
dN/d£ > (mode G in Fig. O?) are unstable (similar to 
the conventional Vakhitov-Kolokolov criterion ^3). In 
Fig. |2f we show the dynamics of an unstable solution (it 
corresponds to mode G in Fig. ^j^) which displays trans- 
formation of an asymmetric mode into a pulsating sym- 
metric distribution at time T w 150. By comparing the 
initial, T = 0, and final, T — 1000, distributions which 
are shown in the insets in Fig.|3f, one can suggest that the 
antisymmetric mode was transformed into a symmetric 
state which is localized about the local minimum of the 
potential i^(X) and nonlinearity G{X). In order to check 
this conjecture we estimated the energy of the emerging 
mode (it is £ — —7.66) and compared its shape with 
the shape of the stationary mode obtained for the same 
energy. The result depicted by the dashed line in the 
right inset of Fig. [Sf, illustrating close similarity between 
slightly oscillating dynamical solution and the stationary 
localized mode, strongly supports the above supposition. 



IV. CONCLUSION 

We have shown that bosonic component of quasi-one- 
dimensional boson-fermion mixtures, loaded in relatively 
deep optical latices can be described by the NLS equa- 
tion with a periodic potential and a periodic nonlinear- 
ity. The fermonic component modifies the effective lat- 
tice for bosons and originates modulation of the interac- 
tion among bosons. In such a system there exist local- 
ized mode solutions, with properties very different form 
those of the known models with homogeneous nonlinear- 
ity and/or potential. In particular, we established regions 
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of parameters, where existence of the locahzed modes re- 
quires a minimal number of bosons and is not Hmited 
by some upper bound (these phenomena resemble prop- 
erties of the quintic NLS equation ^^). The respective 
dependence of the number of the mode particles on the 
energy, unlike in any other known NLS models with peri- 
odic coefficients, has a zig-zag behavior, originated by the 
set of bifurcations. We have also found situations where 
no small amplitude gap solitons can exist near either of 
gap edges, thus making the both Bloch states bordering 
the gap to be modulationally stable, and where not for 
all energies in a gap solitary wave solutions are available. 
Most of the symmetric modes were found to be dynami- 
cally stable. 
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while panel f shows the dynamics of mode G (in the insets 
the initial, T — 0, and final, T = 1000, shapes are shown by 
solid lines; the numerically calculated shape of the mode with 
£ = —7.66 is depicted in the right-hand inset in panel (f) by 
the dashed curve). 
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